/*
===========================================================================
This source file is part of DyLab (Dynamic Laboratory)
For the latest info, see http://dylab-modules.googlecode.com

Copyright (c) 2006-2008 Lukas Krejci
(krejci.lukas@volny.cz)

This file is part of DyLab.

    DyLab is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    DyLab is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with DyLab.  If not, see <http://www.gnu.org/licenses/>.
===========================================================================
*/

//===========================================================================================
#include <dylab/utils/property/PropertyTree.h>
#include <dylab/utils/types/dynamicCast.h>
#include <dylab/utils/topology/TetraComplex.h>
#include <dylab/utils/math/Math.h>
#include <dylab/utils/log/Log.h>
#include "dotT4Obj_ModelStructureLoader.h"

//===========================================================================================
const dylab::String dylab::dotT4Obj_ModelStructureLoader::COMPONENT_NAME = "t4obj";

//===========================================================================================
dylab::dotT4Obj_ModelStructureLoader::dotT4Obj_ModelStructureLoader()
: ComModelStructureLoader(COMPONENT_NAME), mBoundaryMeshScale(DYLAB_R(1)), mBoundaryMeshTranslate(Vector3::ZERO), 
mNodeStep(16), mNodeMass(DYLAB_R(0.1)), mDebugMode(false)
{
}
//-----------------------------------------------------------------------------
void dylab::dotT4Obj_ModelStructureLoader::loadFromPropertyTree(const PropertyTree & _pt)
{
    // retrieve the properties
    mBoundaryMeshScale = _pt.getNodeValueT("boundary_mesh_scale", mBoundaryMeshScale);
    mBoundaryMeshTranslate = _pt.getNodeValueT("boundary_mesh_translate", mBoundaryMeshTranslate);
    mNodeStep = _pt.getNodeValueT("node_step", mNodeStep);
    mNodeMass = _pt.getNodeValueT("node_mass", mNodeMass);
    mDebugMode = _pt.getNodeValueT("debug_mode", mDebugMode);

    // build the structure
    _buildModel(_pt.getNodeValue("boundary_mesh"));
}
//-----------------------------------------------------------------------------
void dylab::dotT4Obj_ModelStructureLoader::onCreateVertex(const Vector3 & _v)
{
    if (mDebugMode)
    {
        DYLAB_LOGF("create vertex: %1%", (_v * mBoundaryMeshScale + mBoundaryMeshTranslate));   
        getModelToLoad()->getBoundaryMeshInterface()->createVertex(_v * mBoundaryMeshScale + mBoundaryMeshTranslate, Vector3::ZERO, 0, 1);
        return;
    }

    // create new node
    NodeSimplex * node = new NodeSimplex(mBoundaryFaces.simplicesReader(NodeSimplex::ORDER).getCount(), &mBoundaryFaces);
    DYLAB_ASSERT(node != NULL);

    // set its attributes
    node->position = _v * mBoundaryMeshScale + mBoundaryMeshTranslate;

    DYLAB_LOGF("create vertex: %1%", node->position);   

    // merge into the bounding box
    mBoundaryMeshBoundingBox.merge(node->position);
}
//-----------------------------------------------------------------------------
void dylab::dotT4Obj_ModelStructureLoader::onSetVertexNormal(uint32_t _index, const Vector3 & _n)
{
    DYLAB_LOGF2("set vertex %2% normal: %1%", _n, _index);

    if (mDebugMode)
    {
        getModelToLoad()->getBoundaryMeshInterface()->vertexReader()[_index]->normal() = _n;
        return;
    }

    // set node normal vector
    mBoundaryFaces.getNode(_index)->normal = _n;
}
//-----------------------------------------------------------------------------
void dylab::dotT4Obj_ModelStructureLoader::onSetVertexTexture(uint32_t _index, real_t _u, real_t _v, real_t _w)
{
    DYLAB_LOGF4("set vertex %1% texture: %2% %3% %4%", _index, _u, _v, _w);
    
    if (mDebugMode)
    {
        ComModelStructure::IBoundaryMesh::Vertex * v = getModelToLoad()->getBoundaryMeshInterface()->vertexReader()[_index];
        v->setTextureU(_u);
        v->setTextureV(_v);
        return;
    }

    NodeSimplex * node = mBoundaryFaces.getNode(_index);
    node->u = _u;
    node->v = _v;
}
//-----------------------------------------------------------------------------
void dylab::dotT4Obj_ModelStructureLoader::onCreateFace(uint32_t _v0, uint32_t _v1, uint32_t _v2)
{
    DYLAB_LOGF3("create face: %1% %2% %3%", _v0, _v1, _v2);

    if (mDebugMode)
    {
        getModelToLoad()->getBoundaryMeshInterface()->createFace(_v0, _v1, _v2);
        return;
    }
    
    // craete face
    Face * face = new Face(&mBoundaryFaces, _v0, _v1, _v2);
    DYLAB_ASSERT(face != NULL);    
}
//-----------------------------------------------------------------------------
void dylab::dotT4Obj_ModelStructureLoader::_buildModel(const String & _boundaryMeshFile)
{
    // parse the given boundary mesh file
    WavefrontObjFormat * parser = new WavefrontObjFormat(this);
    DYLAB_ASSERT(parser != NULL);
    parser->parseFile(_boundaryMeshFile);
    delete parser;

    DYLAB_LOGF("Boundary mesh faces: %1%", mBoundaryFaces.simplicesReader(Face::ORDER).getCount());
    DYLAB_LOGF("boundary mesh bb: %1%", mBoundaryMeshBoundingBox);

    // debug mode - no teta generation
    if (mDebugMode)
    {
        DYLAB_LOG("dotT4Obj_ModelStructureLoader: Debug Mode on!");
        return;    
    }

    // retrieve bounding box extents
    Vector3 extents = mBoundaryMeshBoundingBox.getExtents();
    DYLAB_LOGF("boundary mesh extents: %1%", extents);
    uint32_t x = static_cast<uint32_t>(extents.x());
    uint32_t y = static_cast<uint32_t>(extents.y());
    uint32_t z = static_cast<uint32_t>(extents.z());

    // compute tetra-box extents
    uint32_t tbX, tbY, tbZ;
    tbX = x + (mNodeStep - (x % mNodeStep));
    tbY = y + (mNodeStep - (y % mNodeStep));
    tbZ = z + (mNodeStep - (z % mNodeStep));
    DYLAB_LOGF3("tetra-box extents: %1% %2% %3%", tbX, tbY, tbZ);

    // generate tetra-box
    _generateBox(mBoundaryMeshBoundingBox.getMinPoint(), tbX, tbY, tbZ);

    // cull-out all degenerated tetrahedra
    //_removeDegeneratedTetrahedra();

    // first stage: cull all tetrahedra that lies completely outside the boundary
    _cullOuterTetrahedra(CM_TEST_BOUNDARYFACES);

    // cull-out all degenerated tetrahedra
    _removeDegeneratedTetrahedra();

    // second stage: merge all boundary vertices into tetrahedra nodes, effectively subdividing boundary tetrahedra
    _subdivideTetrahedraByBoundaryVertices();

    // third stage: subdivide tetrahedra by intersection of its edge and boundary face 
    _subdivideTetrahedraByBoundaryFaces();

    // fourth state: subdivide tetrahedra by intersection of its face and boundary edge
    _subdivideTetrahedraByBoundaryEdges();

    // fifth stage: cull all tetrahedra, whos centers lies outside the boundary
    _cullOuterTetrahedra(CM_TEST_BARYCENTER);

    // and finally, generate boundary mesh
    ComModelStructure::IBoundaryMesh * boundaryMesh = getModelToLoad()->getBoundaryMeshInterface();
    boundaryMesh->generateFromTetrahedralVolume();

    // TODO: set (generate?) the uv-texture mapping for each boundary vertex
    for (ComModelStructure::IBoundaryMesh::VertexArray::Iterator vertex(boundaryMesh->vertexIterator()); vertex.isValid(); vertex++)
    {
        vertex->setTextureU(0);
        vertex->setTextureV(1);
    }
}
//-----------------------------------------------------------------------------
void dylab::dotT4Obj_ModelStructureLoader::_generateBox(const Vector3 & _start, uint32_t _width, uint32_t _height, uint32_t _depth)
{
    int32_t nodeCount, cubeCount, i;

    // compute number of nodes in each direction
    int32_t xNodes = _width / mNodeStep + 1;
    int32_t yNodes = _height / mNodeStep + 1;
    int32_t zNodes = _depth / mNodeStep + 1;

    // store the number of all box nodes and cube count (each cube is subdivided into 5 tetrahedrons)
    nodeCount = xNodes * yNodes * zNodes;
    cubeCount = (xNodes - 1) * (yNodes - 1) * (zNodes - 1);

    // retrieve mass points interface
    ComModelStructure::IMassPoints * massPoints = getModelToLoad()->getMassPointsInterface();
    ComModelStructure::IMassPoints::MassPoint * mp;

    //
    ComModelStructure::IMassPoints::MassPoint::id_t firstMP = massPoints->createMassPoints(nodeCount);

    // first pass is going through all box nodes initializing its attributes based on supplied parameters
    for (int32_t y = 0; y < yNodes; y++) 
        for (int32_t z = 0; z < zNodes; z++) 
            for (int32_t x = 0; x < xNodes; x++) 
    {
        // compute node index based on its position within the box
        i = x + z * xNodes + y * (zNodes * xNodes);

        // create new mass point
        mp = massPoints->massPointReader().get(firstMP + i);

        // set mass point position in the box
        Vector3 pos = Vector3(x, y, z).mul(mNodeStep) + _start;
        mp->currState().position = pos;

        // set mass point texture coordinates based on mass point position in terms of the box - TODO: implement better texture mapping
        // TODO
        //mp->setU(float(x + z * xNodes) / float(xNodes * zNodes)).setV(float(y) / float(yNodes - 1));

        // set mass point mass, all nodes have the same mass
        mp->setMass(mNodeMass);

        // set boundary flag
        if ((x == 0) || (x == xNodes - 1) || (y == 0) || (y == yNodes - 1) || (z == 0) || (z == zNodes - 1))
            mp->flags().set(ComModelStructure::IMassPoints::MassPoint::F_BOUNDARY);
    }

    // utility variables
    int32_t n, firstNode, xzNodes;
    xzNodes = xNodes * zNodes;
    bool rotateCube = true;

    // retrieve tetrahedral volume interface
    ComModelStructure::ITetrahedralVolume * tetraVol = getModelToLoad()->getTetrahedralVolumeInterface();

    // second pass goes through indiviadual cubes creating 5 tetrahedrons for each one of them
    for (i = 0, n = 1; i < cubeCount; i++, n++)
    {
        if ((n % xNodes) == 0) { // last node of the row => new row (z++)
            n++;
            if ((xNodes % 2) == 1) // odd number of xNodes => extra switch of the rotateCube flag when entering new row
                rotateCube = !rotateCube;
        }
        if (((n + xNodes - 1) % xzNodes) == 0) // first node on the last row => new level (y++)
            n += xNodes;

        // index of the cube first node (upper-left node at the bottom plane)
        firstNode = n - 1;

        if (rotateCube) // this tetrahedron box will be rotated around axis y about 90 degrees
            tetraVol->generateCubeSection(firstNode, firstNode + 1, firstNode + 1 + xNodes, firstNode + xNodes, 
                            xzNodes + firstNode, xzNodes + firstNode + 1, xzNodes + firstNode + 1 + xNodes, xzNodes + firstNode + xNodes);
        else
            tetraVol->generateCubeSection(firstNode + 1, firstNode + 1 + xNodes, firstNode + xNodes, firstNode,
                            xzNodes + firstNode + 1, xzNodes + firstNode + 1 + xNodes, xzNodes + firstNode + xNodes, xzNodes + firstNode);
 
        // swap the the cube rotation flag
        rotateCube = !rotateCube;
    }
}
//-----------------------------------------------------------------------------
void dylab::dotT4Obj_ModelStructureLoader::_cullOuterTetrahedra(cullmethod_t _cullMethod)
{
    // retrieve mass points & tetrahedral volume interface
    ComModelStructure::IMassPoints * massPoints = getModelToLoad()->getMassPointsInterface();
    ComModelStructure::ITetrahedralVolume * tetraVol = getModelToLoad()->getTetrahedralVolumeInterface();

    // for cull info logging
    uint32_t culledTotal = 0;

    // walk through all tetrahedra
    for (ComModelStructure::ITetrahedralVolume::TetrahedronMap::Iterator tetra(tetraVol->tetrahedronIterator()); tetra.isValid(); )
    {        
        // flag determining, whether the tetrahedron is outer one
        bool outer = true;

        switch (_cullMethod)
        {
        case CM_TEST_APEXES:
            // test all tetrahedron nodes, if they lie inside the boundary        
            for (TetrahedronGeometry::VertexIdxIterator vertexIdx; vertexIdx.isValid(); vertexIdx++)
            {
                // the node is inside
                if (!_isOutsideBoundary(tetra->getMassPoint(vertexIdx)->currState().position))
                {   // so the tetrahedron is not outer
                    outer = false;
                    break;
                }
            }
            break;
        case CM_TEST_BOUNDARYFACES:
            // test all tetrahedron nodes, if they lie inside the boundary        
            for (TetrahedronGeometry::VertexIdxIterator vertexIdx; vertexIdx.isValid(); vertexIdx++)
            {
                // the node is inside
                if (!_isOutsideBoundary(tetra->getMassPoint(vertexIdx)->currState().position))
                {   // so the tetrahedron is not outer
                    outer = false;
                    break;
                }
            }
            if (outer)
            { 
                // retrieve tetrahedron vertices
                TetrahedronGeometry tetraGeom = tetra->getGeometry();

                // walk through all boundary faces
                for (TriComplexT::Simplices::Iterator simplex(mBoundaryFaces.simplicesIterator(Face::ORDER)); simplex.isValid(); simplex++)
                {
                    // retrieve face vertices
                    Face * face = dynamicCast<Face>(simplex());
                    TriangleGeometry3 faceGeom(mBoundaryFaces.getNode(face->getNodeID(TriangleGeometry3::V_0))->position, 
                                               mBoundaryFaces.getNode(face->getNodeID(TriangleGeometry3::V_1))->position,
                                               mBoundaryFaces.getNode(face->getNodeID(TriangleGeometry3::V_2))->position);

                    // find the tetrahedron, it is located in
                    for (TetrahedronGeometry::EdgeIdxIterator edgeIdx; edgeIdx.isValid(); edgeIdx++)
                    {
                        // test, if the the face intersects the edge  
                        if (faceGeom.computeEdgeIntersection(tetraGeom.getEdgeGeometry(edgeIdx)) == Math::INTERSECTION_POINT)
                        {
                            // the tetahedron intersects with the triangle -> the tetrahedron is not outer
                            outer = false;
                            break;
                        }
                    }        

                    if (!outer)
                        break;
                }
            }
            break;
        case CM_TEST_BARYCENTER:
            {
            // compute tetrahedron's center
            Vector3 center = tetra->getGeometry().computeCenter();             

            // debug verifying
            if (tetra->getGeometry().classifyPoint(center).first != Math::POINTCLS_INNER) 
                DYLAB_LOGF4("!!!Tetrahedron (%1%) barycenter (%2%) not in (degenerated=%3%, value=%4%)!", 
					tetra->getGeometry(), center, (uint32_t)tetra->getGeometry().isDegenerated(), tetra->getGeometry().computeValue());

            // test the berycenter againts the boundary
            outer = _isOutsideBoundary(center);            
            }
            break;
        }

        // if the tetrahedron is completely outside the boundary
        if (outer) // then remove it
        {
            ComModelStructure::ITetrahedralVolume::Tetrahedron::id_t tetraID = ~tetra;
            tetra++;
            tetraVol->destroyTetrahedron(tetraID);
            culledTotal++;
        }
        else
        {
            tetra++;
        }
    }

    // log cull info
    DYLAB_LOGF("Culled tetrahedra total: %1%", culledTotal);    
}
//-----------------------------------------------------------------------------
bool dylab::dotT4Obj_ModelStructureLoader::_isOutsideBoundary(const Vector3 & _p)
{
    // construct arbitrary ray from the node position
    Vector3 rand = Vector3::random();
    Ray3 ray(_p, rand - _p);
    ray.direction().normalise();

    // test the vector for intersection with all boundary face
    bool outside = true;

    //
    Vector3 currIP, lastIP;
    bool    lastIPPresent = false;    

    for (TriComplexT::Simplices::Iterator simplex(mBoundaryFaces.simplicesIterator(Face::ORDER)); simplex.isValid(); simplex++)
    {
        Face * face = dynamicCast<Face>(simplex());
        TriangleGeometry3 faceGeom(mBoundaryFaces.getNode(face->getNodeID(TriangleGeometry3::V_0))->position,
                                   mBoundaryFaces.getNode(face->getNodeID(TriangleGeometry3::V_1))->position,
                                   mBoundaryFaces.getNode(face->getNodeID(TriangleGeometry3::V_2))->position);        

        if (faceGeom.computeRayIntersection(ray, currIP) == Math::INTERSECTION_POINT)
        {
            outside = !outside; // even number of intersection -> outside, odd -> inside                

            if (lastIPPresent && currIP.almostEqual(lastIP))            
                outside = !outside;
                        
            lastIP = currIP;
            lastIPPresent = true;
        }
    }

    // return the point classification result flag
    return outside;
}
//-----------------------------------------------------------------------------
void dylab::dotT4Obj_ModelStructureLoader::_removeDegeneratedTetrahedra()
{
    // retrieve mass points & tetrahedral volume interface
    ComModelStructure::IMassPoints * massPoints = getModelToLoad()->getMassPointsInterface();
    ComModelStructure::ITetrahedralVolume * tetraVol = getModelToLoad()->getTetrahedralVolumeInterface();

    // for info logging
    uint32_t degenerated = 0;

    // walk through all tetrahedra
    for (ComModelStructure::ITetrahedralVolume::TetrahedronMap::Iterator tetra(tetraVol->tetrahedronIterator()); tetra.isValid(); )
    {        
        // if the tetrahedron is degenerated
        if (tetra->getGeometry().isDegenerated()) // then remove it
        {
            DYLAB_LOGF("$$$ Degenerated tetrahedron found: %1%, removing it.", tetra->getGeometry());
            ComModelStructure::ITetrahedralVolume::Tetrahedron::id_t tetraID = ~tetra;
            tetra++;
            tetraVol->destroyTetrahedron(tetraID);
            degenerated++;
        }
        else
        {
            tetra++;
        }
    }

    // log info
    DYLAB_LOGF("Degenerated removed: %1%", degenerated);    
}
//-----------------------------------------------------------------------------
dylab::uint32_t dylab::dotT4Obj_ModelStructureLoader::_subdivideTetrahedraByBoundaryVertices()
{
    // retrieve mass points & tetrahedral volume interface
    ComModelStructure::IMassPoints * massPoints = getModelToLoad()->getMassPointsInterface();
    ComModelStructure::ITetrahedralVolume * tetraVol = getModelToLoad()->getTetrahedralVolumeInterface();

    // for log info
    uint32_t subdividedTetraCount = 0;

    // walk through all boundary vertices
    for (TriComplexT::Simplices::Iterator simplex(mBoundaryFaces.simplicesIterator(NodeSimplex::ORDER)); simplex.isValid(); simplex++)
    {
        Vector3 p = dynamicCast<NodeSimplex>(simplex())->position;
        bool found = false;

        // find the tetrahedron, it is located in
        for (ComModelStructure::ITetrahedralVolume::TetrahedronMap::Iterator tetra(tetraVol->tetrahedronIterator()); tetra.isValid(); )
        {
            // test, if the subdivision point is located insinde the tetrahedron
            TetrahedronGeometry::pointcls_t pointCls = tetra->getGeometry().classifyPoint(p);
            if ((pointCls.first != Math::POINTCLS_OUTER) && (pointCls.first != Math::POINTCLS_VERTEX))
            {
                // tetrahedron has been found, subdivide it and break the loop
                DYLAB_LOGF3("Vertex %1% penetrates tetrahedron %2% at %3%", p, tetra->getGeometry(), pointCls);
                found = true;                

                switch (pointCls.first)
                {
                case Math::POINTCLS_INNER:
                    {
                        ComModelStructure::ITetrahedralVolume::Tetrahedron::id_t tetraID = ~tetra; 
                        tetra++;
                        tetraVol->subdivideTetrahedron(tetraID, p); 
                    }
                    break;
                case Math::POINTCLS_EDGE:
                    {
                        ComModelStructure::ITetrahedralVolume::Edge::id_t edgeID = tetra->getEdgeID(pointCls.second.edgeIdx);
                        tetra++;
                        tetraVol->subdivideEdge(edgeID, p);
                    }
                    break;
                case Math::POINTCLS_FACE:
                    {
                        ComModelStructure::ITetrahedralVolume::Triangle::id_t triID = tetra->getTriangleID(pointCls.second.faceIdx);
                        tetra++;
                        tetraVol->subdivideTriangle(triID, p);
                    }
                    break;
                }

                _removeDegeneratedTetrahedra();                
                subdividedTetraCount++;
                break;
            }            
            else
            {
                tetra++;
            }
        }

        if (!found)
            DYLAB_LOGF("*** Vertex %1% penetrates no tetrahedron!", p);
    }

    // log subdivision info
    DYLAB_LOGF("Subdivided tetrahedra total: %1%", subdividedTetraCount);

    return subdividedTetraCount;
}
//-----------------------------------------------------------------------------
dylab::uint32_t dylab::dotT4Obj_ModelStructureLoader::_subdivideTetrahedraByBoundaryFaces()
{
    // retrieve mass points & tetrahedral volume interface
    ComModelStructure::IMassPoints * massPoints = getModelToLoad()->getMassPointsInterface();
    ComModelStructure::ITetrahedralVolume * tetraVol = getModelToLoad()->getTetrahedralVolumeInterface();

    // for log info
    uint32_t subdividedEdgeCount = 0;

    // walk through all boundary faces
    for (TriComplexT::Simplices::Iterator simplex(mBoundaryFaces.simplicesIterator(Face::ORDER)); simplex.isValid(); simplex++)
    {
        // retrieve face vertices
        Face * face = dynamicCast<Face>(simplex());
        TriangleGeometry3 faceGeom(mBoundaryFaces.getNode(face->getNodeID(TriangleGeometry3::V_0))->position, 
                                mBoundaryFaces.getNode(face->getNodeID(TriangleGeometry3::V_1))->position,
                                mBoundaryFaces.getNode(face->getNodeID(TriangleGeometry3::V_2))->position);

        // gather all current edges
        ComModelStructure::ITetrahedralVolume::EdgeList edgeList(false);
        for (ComModelStructure::ITetrahedralVolume::EdgeMap::Iterator edge(tetraVol->edgeIterator()); edge.isValid(); edge++)
            edgeList.insertAtEnd(edge());

        // walk through all edges gathered for this face
        for (ComModelStructure::ITetrahedralVolume::EdgeList::Iterator edge(edgeList); edge.isValid(); edge++)
        {
            // intersection point
            Vector3 ip;

            // retrieve geoemtry of the edge
            EdgeGeometry3 edgeGeom = edge->getGeometry();

            // test, if the the face intersects the edge           
            if (faceGeom.computeEdgeIntersection(edgeGeom, ip) == Math::INTERSECTION_POINT)
            {
                // classify the intersection point with respect to the face
                EdgeGeometry3::pointcls_t pointCls = edgeGeom.classifyPoint(ip, DYLAB_R(100) * Math::EQUALITY_TOLERANCE, DYLAB_R(1));

                if (pointCls.first == Math::POINTCLS_INNER)
                {                
                    DYLAB_LOGF3("Vertex %1% penetrates edge %2% at %3%", ip, edgeGeom, pointCls);

                    // subdivide the edge
                    tetraVol->subdivideEdge(edge->getId(), ip); 
                    subdividedEdgeCount++;

                    // remove degenerated tetrahedra
                    _removeDegeneratedTetrahedra();                
                }            
            }
        }        
    } 

    // log subdivision info
    DYLAB_LOGF("Subdivided edges total: %1%", subdividedEdgeCount);

    return subdividedEdgeCount;
}
//-----------------------------------------------------------------------------
dylab::uint32_t dylab::dotT4Obj_ModelStructureLoader::_subdivideTetrahedraByBoundaryEdges()
{
    // retrieve mass points & tetrahedral volume interface
    ComModelStructure::IMassPoints * massPoints = getModelToLoad()->getMassPointsInterface();
    ComModelStructure::ITetrahedralVolume * tetraVol = getModelToLoad()->getTetrahedralVolumeInterface();

    // for log info
    uint32_t subdividedTriangleCount = 0;
    uint32_t pass = 0;
    bool restart;

    do
    {
        uint32_t passSub = 0;
        pass++;
        restart = false;
    // walk through all boundary faces
    for (TriComplexT::Simplices::Iterator simplex(mBoundaryFaces.simplicesIterator(EdgeSimplex::ORDER)); simplex.isValid(); simplex++)
    {
        // retrieve face vertices
        EdgeSimplex * edge = dynamicCast<EdgeSimplex>(simplex());
        EdgeGeometry3 edgeGeom(mBoundaryFaces.getNode(edge->getNodeID(EdgeGeometry3::V_0))->position, 
                               mBoundaryFaces.getNode(edge->getNodeID(EdgeGeometry3::V_1))->position);

        // gather all current edges
        ComModelStructure::ITetrahedralVolume::TriangleList triangleList(false);
        for (ComModelStructure::ITetrahedralVolume::TriangleMap::Iterator triangle(tetraVol->triangleIterator()); triangle.isValid(); triangle++)
            triangleList.insertAtEnd(triangle());

        // find the tetrahedron, it is located in
        for (ComModelStructure::ITetrahedralVolume::TriangleList::Iterator triangle(triangleList); triangle.isValid(); triangle++)
        {
            // intersection point
            Vector3 ip;
            TriangleGeometry3 faceGeom = triangle->getGeometry();

            // test, if the the face intersects the triangle
            if (faceGeom.computeEdgeIntersection(edgeGeom, ip) == Math::INTERSECTION_POINT)
            {
                // classify the intersection point with respect to the face
                TriangleGeometry3::pointcls_t pointCls = faceGeom.classifyPoint(ip);

                if ((pointCls.first != Math::POINTCLS_OUTER) && (pointCls.first != Math::POINTCLS_VERTEX))
                {                
                    DYLAB_LOGF3("Vertex %1% penetrates triangle %2% at %3%", ip, faceGeom, pointCls);

                    switch (pointCls.first)
                    {
                    case Math::POINTCLS_INNER:                        
                        tetraVol->subdivideTriangle(triangle->getId(), ip);                         
                        break;
                    case Math::POINTCLS_EDGE:                        
                        tetraVol->subdivideEdge(triangle->getEdgeID(pointCls.second.edgeIdx), ip);                        
                        break;
                    }

                    _removeDegeneratedTetrahedra();                
                    subdividedTriangleCount++;
                    passSub++;
                    restart = true;
                    break;
                }            
            }
        }    
    }
    DYLAB_LOGF2("^^^ pass %1% subdivided: %2%", pass, passSub);
    } while (restart);

    // log subdivision info
    DYLAB_LOGF("Subdivided triangles total: %1%", subdividedTriangleCount);

    return subdividedTriangleCount;
}
//-----------------------------------------------------------------------------